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Instead, we have uncovered a set of grid cells that 
represent a rapid change in shape. When we encounter a 
grid cell of .this type, then we deform the shape in the 
grid cell to conform to any adjacent grid cell and assign 
5 it to the equivalence class that forces the least amount 
of deformation. 

Two equivalence class shapes may fail to connect in 
a region where they should. When this happens, we deform . 
the respective shape in those boundary grid cells until 
10 they match the tri- linear interpolation function in each 
grid cell. 

We now discuss the interface between two grid cell 
equivalence classes . If a grid cell face separates the 
two equivalence classes, then we call the interface 

15 "sharp", otherwise the interface is "diffuse". If the 
interface is "sharp" , then there exists a ID boundary 
separating the two shapes. We compute this boundary by 
equating the. two quadratic surface expressions and 
solving. This is not always smooth, so it will be 

20 necessary to blend a narrow region on* each side of the 

intersection curve.. If the interface is diffuse, then we 
have a 2D region of overlap between the two shapes in 
which the smoothing can be naturally accommodated. We fit 
a shape to the overlap. We close this topic by checking 

25 this analysis against two adjacent grid cells that have 
different conic section fibre bundle shape 
approximations. We conclude that the analysis is valid, 
since tri-linear interpolation is continuously 
differentiate on the grid cell face that separates the 

30 two adjacent grid cells. 

ISfow we overlay the patchwork shape on top of the 
• original implicit surface S and look for chances to 
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reduce the misfit. We will do this by approximating the 
volume as a 3D Riemann integral. Our misfit reduction 
procedure is as follows. 

Misfit reduction algorithm 

1. We project the noisy implicit surface onto 
the patchwork smooth surface. We identify 
connected regions of the implicit surface 
that are responsible for significant misfit 
and try to reduce these regions first. 

2 . We anticipate 'multiple misfit regions 
associated to each smooth shape. The most 
complicated element in a misfit description 
is the representation of the outermost ID 
boundary of the misfit. We use the outermost 
ID boundary to define a Riemann sum 
approximation to the misfit. Consequently, .we 
will ignore (initially, at least) low 
amplitude but- high frequency nearby features. 

.3. Define .height field (relative to the smooth 
approximant) contours at each critical point 
in the noisy implicit surface. Project these 
contours onto the smooth approximant. 

4. Smooth the set of projected contours. For 
each contour in the parameterization we 
identify the critical points on that contour. 
We triangulate the annular region bounded 
between consecutive contours . 

5. Now we fit a f rustum-like collar to the base 
of the region where the misfit attaches to 
the background. 

6. * The cumulative shape resembles a telescope, 

hence the name "telescopic deformation". The 
process accommodates tendril-like shapes, in 
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the sense that surfaces can recursively 
support branches off other earlier branches. 

Figures 28a-c illustrate an example of the misfit 
5 reduction process . Suppose that we want to reduce the 
misfit in approximating a tendril shape by a plane, see 
Figure 28a. We subdivide the tendril using height field 
planes that approximately set to the tendril's critical 
points . 

10 We have generated a sequence of truncated collars 

that reduce the misfit. We are not done, however, since 
this approximation is not dif f erentiable at any point on 
. a collar ring separating two truncated cones . We remedy 
this deficiency by replacing a thin swath about the ring 
15 with a thin swath about the equator of a sphere. 

Figures 29a-c illustrates an example of blending a 
non-diff erentiable join of two collars. The idea is to 
. map given instance (Figure 29a) to the one situation that 
we understand how to fix (Figure 29b) . -The one situation 
20 that we understand hdw to fix is an isosceles right 

triangle with an inscribed circle such that the center of 
the circle is the centroid of the triangle. We apply the 
linear transformation- shown in Figure 29c to the 
isosceles right triangle in panel shown in Figure 29b to 
25 the exterior and interior triangles - corresponding to 0i 
and 6 2 , respectively - in Figure 29a. 

The misfit reduction procedure is as follows. 

. l.We construct height field contours on the noisy- • 
30 surface S. Remember that the height field is taken 

• normal to the smooth shape approximation S* . 
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2. Working top down, we project each contour to S*. It 
is possible that the projection of successive 
contours cross . - 

3. We approximate each contour projection in two" steps. 
5 First we locate all ID critical points - minima, 

maxima, and inflection points. Second we fit a 

smooth approximation to each section of the 

projected contour; where a section is delimited by 

critical points. 
10 4 . We triangulate the interior of each pro j ected 

contour, and then construct tetrahedra using 

consecutive contours . 
5 . We approximate the cap/cup cbrresponding to the 

innermost ring of the bull's eye with the ring's 
15 closure in the plane. 

We next consider how to describe a branching 
structure of hills that separate a region into multiple 
valleys. Figure 30 shows as a geological example a water 
20 breach as indicated by the white arrow. 

Here is a compact description of this branching 
structure. 

1. Find the ridge edge in the branching structure. The 
25 ridge edge is a connected (therefore degenerate) set 

of critical points . 

2. We assume that as typical cross section is a conic 
section, e.g., a parabola, an ellipse, etc. 

3 . We allow for stretching and contracting of the cross 
30 - section. We call a region over which the typical 

cross section stretches or contracts a "transition 
zone".' We provide linear interpolation between the 
start and stop position cross sections. 
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4. We anticipate a valley that runs perpendicularly to 
the ridge edge. For the case of a branching 
structure, we treat a valley of this type as a 
misfit. The following NASA aerial reconnaissance 
5 image contains active breaches in the ridge 

structure..' 

We consider a more challenging example. Figure 31 is 
a NASA Shuttle Mission photograph of the Richat Structure 

10 in Mauritania. NASA believes that most likely explanation 
of the origin of this structure is that it is uplifted 
rock sculpted by erosion. NASA also says that there is no 
widely accepted theory that explains why the Richat 
Structure is nearly circular. 

15 The dominant shape in the Richat Structure is a 

nested sequence of toroidal structures. Because of 
erosion, some of the toroidal structures exist as 
sections of a torus only. Ecker argues on pg. 4-5 of his 
lecture notes that a torus .satisfies; mean curvature flow 

20 with respect to its mean curvature. This flow is 

singular, in the sense that in finite time the circular 
cross section shrinks to a point and hence the torus 
shrinks to a circle. This is a second example of the fact 
that mcf evolution need not preserve the dimension of a 

25 shape. (The other example is an uncapped cylinder whose 
radius r(t) is given by the formula 

r(jt) = y[r 2 (P)-2t 7 provided that *e (Q,— ^-). In finite time the 

uncapped cylinder collapses to its axis of symmetry. 
With respect to known techniques, Rouby et al. 
30 describe a parameterization algorithm for triangulated 
surfaces. See, D. Rouby, H. Xiao, J. Suppe, 3D 
Restoration of Completely Folded and Faulted Surfaces 
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Using Multiple' Unfolding Mechanisms , AAPG Bulletin 84(6), 
June 2000, pg. 805-829, This method involves projecting 
the ■ triangulation onto the (x,y) co-ordinate plane. It is 
not clear how to compensate for triangles that have a 
5 near singular projection. Also, the projection of 

multiple triangles • onto a common plane can lead to a 
complicated overlapping and slivers. These are the sort 
of numerical instabilities- that we want to avoid and 
'which this algorithm does avoid. 

10 

A description will now be give for 3D conformal grid 
generation and reconstruction of shape from disconnected 
pieces, according to the invention. 

Shape identification improves the performance of our 
15 incremental 3D conformal grid generator. Given .a sub- 
volume V with boundary 9v = S, we uniformly approximate S 
as a collection of trimmed ideal elementary shapes {Em} . 
Then we enclose S in a tight sequence of rectangular 
prisms { P*} . such that the "vertical" edges of the prism 
. 2 0 are normal to the ideal elementary shape approximants . 

We generate a 3D conformal grid in 3 parts '. 

1. 9P k n S, which is the part between the prism faces 
25 and the measured surface S. 

2 . The bounding box B ro of each trimmed ideal elementary 

shape E TO . 

3. A prism P* £ V such that P* n B m = 0. 

30 We use a proportional spacing correlation scheme to 

grid part #1. In part #2 we use Em's parametric 
expression for normals to grid Bm. In part #3 we grid P* 
in the obvious way. The incremental nature of the 
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procedure follows from the subdivision of the grid into 
regions associated with the prismatic enclosure of S and 
the uniform approximation by trimmed ideal elementary 
shapes . 

5 Figure 32 shows a conformal grid induced from a 

proportionally spaced correlation scheme. We use a 
proportional spacing correlation scheme to define the 
grid iso-surf aces . Here is an ideal application of this 
mechanism (This is an idealization of a cross section 

10' through the Gullfaks field.). We see that the sub-blocks 
have approximately the same z-direction thickness. It is 
irrelevant how the horizontal extents correspond. In 
other words, we can generate either a regular spacing 
grid or an irregular spacing grid. We remark that a 

15 compromise such as a "tartan grid" (also known as a 

rectilinear grid) is more memory efficient than a regular 
grid but still has a convenient algebraic lookup 
function . 

We observe that a correlation scheme implements a 
20 .cheap form of mean curvature flow and if generated 

independently of a parameterization, then the conformal 
grid is a roundabout way to also generate a 
par ame terization. 

25 We conclude this part of the description with a 

discussion of a way to convert an existing 3D Cartesian 
grid to a conformal grid. Figure 33 illustrates a non- 
conformal 3D Cartesian grid. As is well known, it is non- 
trivial to trim grid cells that cross volume . boundaries . 

30 Here the background grid shown in dashed lines, crosses 
the overburden. It is clear in this diagram that the, 
background grid cells do not conform to lithological 
boundaries . 
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We remedy these deficiencies in the following way. 

1. We attach the volume of interest to a spatial 
frame. 

5 2 . We conf ormally grid the volume that lay between 

the overburden and the reservoir. 
3 - For best results we want a transition zone from a 
Cartesian grid to the conformal grid, preventing 
reverberation in the grid across lithological 
10 boundaries . 

Forensic reconstruction, will now be discusses. 
Suppose that an application identifies a set of 
disconnected 2D. patches and would like to construct a 

15 surface from the patch set. We show how mean curvature 
flow enables us to construct such a surface - hence the 
sub- title for this discussion. This technique will be 
seen to be similar in execution to grid generation, which 
is why we discuss it here. 

20 Suppose that we have a collection of 2D patches in 

R 3 . We assume that all of the patches are oriented in a 
consistent manner, which must be the case if all of the 
patches are taken from a common surface. We construct the 
solution using patches or parts of patches that have a 

25 common Gaussian curvature signum. We do not insist that 
the Gaussian curvature signum be constant across the 
patch point set. However, if "it is not constant, then we 
use only V those patches and parts of patches that have the 
same signum to construct a surface. This procedure can 

30 create holes ih the surface. We fill in these holes with 
the remainder of patches that were used only partially. 
Some patches with the opposite signum may not fill in 
. holes, but instead form folds in the surface. We 
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construct a second surface from these patch point sets 
and union the two parts together. 

We remark that the following algorithm is applicable 
to the task of sewing together parts of the bounding 
5 surface of a geobody, e.g., a salt body. 

Our method is as follows. 

1. We assume that the velocity of the flow equals l,so 
10 that a time interval is equal to a distance. 

2. For each patch or region of a patch R that has ( + ) 

signum Gaussian curvature, we define a uniformly 
expanding mean curvature flow with initial 
condition equal to R. 
15 3. These patches are taken from a common surface, so we 

sample each flow at the same time step increments . 

4 . We stop a particular flow when it intersects another 

flow pattern. 

5. Eventually each flow intersects a different flow. 
20 6. It may be that the confluence of two flows is not 

dif ferentiable. If this happens, then we define a 
flow to smooth the sharp angle. Alternatively, we 
can smooth the sharp angle using the Squeeze 
operator defined in the Misfit Reduction 
25 algorithm. 

If a formation contains a tunnel, then we may wish 
to reverse the dissolution process, i.e., fill in the 
tunnel. In this situation we will approximate the tunnel 

30 as a collection of uncapped cylindrical and toroidal 

sections. Ecker's lecture notes show that mean curvature 
flow applied to an uncapped" cylinder will shrink the -3D " 
surface onto its ID axis of symmetry in finite time. 
Smocyzk ,has obtained a comparable result for a torus. His 

35 result shows that ei torus collapses under mean curvature 
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flow to its. inner parabolic circle in finite time. We 
conclude that when applied to a tunnel, mean curvature 
flow will collapse the tunnel to a curve in finite time. 

Figure 34 is an image of the Devil's Potholes, South 
Africa. We can identify two different kinds of pothole 
formations in this image. The circled pothole on the 
right is well described as a cylinder. We can collapse it 
to its axis of symmetry - which is a line - under mean 
curvature flow. The pothole on the left is more complex. 
The top of this region is a curved throat, which is 
followed by a sharp edged cone. Finally there is a 
toroidal section that is turning away from the viewer. We 
know that each section will collapse to a singular edge 
in a finite amount of time. Taking the longest time 
interval needed to collapse a section of this structure, 
we guarantee that the entire pothole will collapse to a 
polygonal path. 

We now discuss a related problem. Given an image of 
a layered sequence, suppose that an unconformity. The 
0 unconformity will erase everything in the image that is 
above the unconformity. How good a reconstruction of the 
eroded interface can I achieve, given just the migrated 
image. As an example, we return to the diagram of a cross 
section of a shallow sequence in Turkmanistan shown in 
5 Figure 3 . We aire interested in the guess of what the 
missing section of the layer looks like. 

Our algorithm assumes that the formation satisfies 
. the following assumptions. 

0 1. There exists a reference layer such that the visible 
part of the sequence is well approximated by a 
single reference surface adaptively sampled 



68 



WO 2004/057375 




PCT/GB2003/005395 



• correlation scheme. In the above schematic, the 
vertically striped layer R is the reference layer. 

2 . We reconstruct the eroded region of a layer by 
evolving the reference interface under mean 
curvature flow until the front joins the 
unconformity. In the above schematic, the 
reconstructed interfaces are shown' as dotted 
outlines above the unconformity . 

While it may not be obvious from Figure 3, conformal 
grid generation is a very short time solution to mean 
curvature flow. In other words, if erosion had not been 
imposed then we would recognize the conformal grid 
generation in the framework close to the reference 
surface and flattening further, away. As an illustration, 
we consider Figure 35, which is a NASA LANDSAT image of 
the Yukon River delta (NASA Geomorphology plate D-9) . We 
focus on the boundary cracks that subdivide the rounded . 
joint at the left top of the image. These boundaries 
delimit curvature flow regions that collectively form the 
bend in the deltaic mass. It has been proposed to use a 
diffusive process. See, J. Davis, S. Marschner, M. Garr, 
M. Levoy, Filling holes in complex surfaces using 
volumetric diffusion , First International Symposium on 3D 
Data Processing, Visualization, and Transmission, Padua, 
Italy, June 19-21, 2002. This reference does not make • 
use of curvature information in the algorithm used. 
Furthermore, the reference doe's not deal with closing 3D 
voids such as tunnel closure. 
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We have discussed how to construct a patchwork 
covering of an implicit surface from first order Taylor 
polynomials. The error term is quadratic vyith 
coefficients that are the principal curvatures defined at 
5 a point that is somewhere within the circle of 

convergence about the expansion point. The radius of the 
circle is chosen so that the error term is smaller than a 
user-defined tolerance. The error estimate is not sharp, 
because it is not obvious where to evaluate the error 
10 term's principal curvature values, i.e., what is the 

definition of a. Since this error estimate is not sharp, 
we supplement the Taylor representation with a 'fast 

technique to evaluate a. 

From curvature analysis we derive a representation 
15 of a surface as a connected sum of trimmed elementary 

shape. Some examples of elementary shapes are (sections 
of) an ellipsoid, hyperboloid, cyclide, and a prism. An 
elementary shape has algebraic expressions for mean 
curvature and Gaussian curvature. These expressions can 

20 • be substituted into the quadratic formula to generate an 
algebraic expression for principal curvature. 
Unfortunately the resulting expression is often so 
complex algebraically that a symbolic algebra package 
£uch as Mathematica® must be used for any manipulation 

25 other than simple evaluation. This is awkward in an 

Engineering environment, so a low order Taylor series 
approximation is preferable. 

Here is an intersection algorithm for implicit 
surfaces. This algorithm assumes that every surface's 

30 signed distance function has been defined on a shared 
grid. 
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Notation 

S x implicit surface #1 
a x a point .on S x 
x parametric definition for a x 

S 2 implicit surface #2 

G 2 a point on E 2 

x 2 parametric definition for a 2 . 

G a shared 3D grid 

C a grid cell in G such that C n SI O ,S2 * 0 

c, a point in C n S x n S 2> 

a x2 a point in C n S x n S 2 . 

Intersection algorithm 

1. Find all grid cells {C k } that surface S d intersects. 
A grid cell is involved if and only if the signed 
distance function is not constant at all grid cell 
corners . 

2 . In this algorithm we want to work with surface 
patches that project invertibly onto a face of grid 
cell. We subdivide the part of a surface contained, 
in a grid cell into the desired patchwork by 
checking the surface normal field restricted to the 
grid cell. We place ( a point on the surface in a 
patch exactly when the unit normal at the ppint is 
within tolerance of being parallel to all other 

0 surface points that are already assigned to- the 

patch. . 

For each S, patch find all S 2 patches such that their 
^corresponding sdf patterns do not rule out intersection. 
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For clarity we assume that each surface contains exactly 
one such patch. 

LrOop$0: 

5 

3 . Let c 0 be the midpoint of the line segment X 0 

connecting two faces of C such that the signum of 
the sdf for S x on one of the connected faces is 
different from that on the other, 
10 4. Let j = 0. 

£iOop$l : 

5. Project c 5 onto S x and S 2 , respectively. 
15 6. Label the projection s ld and s 2j/ respectively. 

7. Either = s 2i or not. 

8. If s^ = s 2:J , then go to Label$2. 

9. If s xi * s 2j then construct the line segment in C k 
connecting s ld to s 2d . 

2 0 10. Let c j+1 be the midpoint of X 3 . +x 

11. It might be the case that c d « c j+lt We expect to see 
this when a very thin hyperbolic neck is contained 
in the same cell as C k . 

12. If c^c^, then we want to jump away from this 

25 hopeless local minimum. We do this by comparing the 

signum of the sdf for s y and the signum of the sdf 
for s 2j to signum of sdf obtained on the grid cell 
faces joined by the line segment X 0 . Moving in the 
direction that differs from the current values, jump 

30 to the midpoint along A, 0 , 

13. " Increment j and return to Loop$l. ■ 

Label $2: 
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We use a Newton algorithm to find another point of 
intersection' 

Loop$3s 

14. Let (x^y^) be a point of intersection for S 1 and S 2 
in C and suppose that (x^y^*) is another point of 
intersection in C. Expanding each surface in a 1 st 
order Taylor series about (x^y^) , we have 

Si (*n * .yu*) = s i (*i2 .yia ) + V5 i (*« . yn ) ■ (^12-^12*^12-^2*) 

Sa(*n*/y»*) = ^2(^2^12) + VS 2 (^ a , 3> 12 ) • (x 12 - x, 2 *, y xi - y ia *) 

15, Subtracting, we have 2 linear equations in two 
unknowns. Solve for the non- trivial solution 

0 = V (S 2 - S t )(* 12 , y«) • (*ia " *i2*> yn - yn*) 



16. Return to Loop$3 until (x u *, y^*) is acceptably 

close to being another point of intersection in C. 
20 17. We compute the intersection curve between (x 12/ y 12 ) 

and (x^*, y 12 *) by repeating the midpoint subdivision 
algorithm in Loop$2 . 

18. Return to Loop$3 and repeat the 1 st order Taylor 
series algorithm in Loop$2 on both ends (x^'Yia) and 

25 (x ia * / y l2 *) until the curve crosses the faces of C k . 

19. Return to Loop$3and continue processing grid cells 
in {C k } until exhaustion. 
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Visualization of an intersection curve is as follows. 

1. Project the intersection curve onto the internal 
cylindrical parameter space. 

2. Facet the cylinder using a constrained Delaunay 
algorithm. 
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3. Invert the mapping, projecting the triangulation of 
the region enclosing the u flattened" intersection 
curve to the parent implicit surface. 

5 We have described a curvature-adaptive method for 

placing sample points on a surface and generating a 
narrow band octree from them. Now we, discuss how spatial 
frames can be used to implement an efficient update 
mechanism. For details of the reappearance of SIGMA' s 
10 hybrid grid-mesh concept, refer to U.S. Patent 

Application Serial No. 09/163,075, incorporated herein by 
reference . 

o 

Each geological feature has its own octree which 
overlays part of the background octree. An overlay can 
15 cover a strictly smaller region of ' the background, and 

this small region can be sampled differently (adaptively) 
from the background. The sampling can be adapted to the 
complexity of the framework region that the overlay 
encloses. The overlay changes as the framework region 
20 changes. A separate spatial frame is assigned to control 
each octree overlay. Overlays are loaded on demand. 
Within a spatial frame the sampling. can be adjusted in an 
"■ adaptive manner. A frame boundary sample node can be part 
of multiple octree overlay. An entity that is outside a 
25 particular spatial frame learns about the frame's 

enclosed shape by interrogating the enclosing spatial 
frame. Given an octree identifier and a location inside 
or on the boundary of the octree's enclosing frame, the 
local cpt returns the cpt and sdf for the -remote octree 
30 overlay. Adjacent frames might support a different set of 
and number of framework surfaces. 

Given a. region of interest, we associate to- each 
spatial frame the point set that the f.rame contains. We 
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represent the set of spatial frames associated with a 
region of interest as a directed acyclic graph. Two 
parent frames that partially overlap will be represented 
in the graph as nodes that point to a common descendant . • 
5 Thus spatial frames constitute a topology graph. 

As an example, a vsp was acquired over part of the 
Zechstein Salt formation. Figures 36a and 3 6b show, for 
reference, the background Zechstein Salt and the region 
in the Zechstein where the vsp was acquired. Figure 37 
10 show the frame graph that ties the vsp region of interest 
to the Zechstein Salt background. We use the inner frame 
to outer frame relationship to . express "A is a boundary 
of B" and "A is logically part of B" . In GQI terminology 
a. spatial frame as defined herein is equivalent to a 
15 gqi_Frame_t. We have extended this functionality two 
ways. First, we maintain logical containment, i.e., 
feature relationships. Second, we allow mixing space and 
time frames. 

Spatial frames are useful in subdividing a faulted 
20 region to match regions subject to different physics. For 
example, Figures 38a, 38b, 39a and 39b illustrate the 
separation of faulted sediments from unfaulted sediments. 
Figure 38a and 38b show the case of a sequence of 
sediments such that some but not all of the sediments 
25 have been faulted, and a simple normal fault is involved. 
Figures 39a and 39b show the case of .a fault network 
whete some of the faults in the network emanate from 
other faults. Figure 39b shows that the unique final 
configuration of spatial frames that partition, the 
30 volumes shown in Figure 39a. Using mathematical 

' induction on the number of faults in a compact region of 
interest, we conclude that it is immaterial in what orde>r 
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spatial frames are assigned to isolate distinct 
sedimentary regions. In other words, the final 
subdi vision is always the same no matter what the order 
of isolation is. 



5 



We think of mean curvature flow as an evolutionary 



process. For simplicity, we assume that any evolutionary 
process is a structure group of dif f eomorphisms of some 
4D fibre bundle. All information regarding the process at 
a moment in time is encoded in a time frame. A time frame 
10 contains a reference to a region of interest that in turn 
is represented as a 3D fibre bundle. We do not demand 
that an evolutionary process . provide .a physically 
plausible explanation of a formation state, rather it is 
enough if the visual expression appears plausible. 
15 Figure .40 illustrates a time-lapse seismic evolution 

(steam injection front tracking, with the left panel 
"before" and the right panel "after".). We discuss now 
how we use spatial frames to define a topology graph. Let 
S be an implicit surface and let {E k } be a family of 
20 elementary shapes that approximate S. Each E k reduces its 
misfit with S by fitting a disjoint region of S to a 
telescopic extrusion of E k . 

We use a different spatial frame to contain each 
telescopic extrusion and attach the collection of these 
25 frames. We assign each elementary shape to its own 

spatial frame and refer to the misfit reduction through 
attachment of the spatial frame list to the parent 
elementary shape frame. Figure 41 shows the reference 
* structure for this set of relationships. " 



We summarize the contents of the new class instances 
that we have defined in this paper. 



30 
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Spatial frame class 

1. Database identifier. 

2. The parent spatial frame. 

3. The list of children spatial frames. 

4. a static hash table that provides the address of 'a 
spatial frame given a database identifier. 

5. The outermost grid cell boundary, given as a pair of 
synchronized lists. . 

a. The first list is a set of UK base grid cells. 
The second list is in 1-1 correspondence with 
the first list. 

b. The second list is an ordered list of 

{left /right, up/down...} steps that define the 
grid cell face sequence that form the boundary. 

6. The boundary curve contents' of the grid cell 
boundary. 

a. Mandatory is a list, of minima, maxima, and 
inflection points, ordered by arc length 
coordinate from the first critical point in 
this list. 

- b. Optional is a list of conic sections that 
approximate the shape of the boundary 
restricted to each consecutive .pair of critical 
points . 

7. Transition zone, given as a grid cell thickness 
surrounding the interior of the outermost boundary.. 
Zero thickness is acceptable. 

8. The list of data base identifiers used to define the 
surfaces that intersect this spatial frame. 

9. The data base identifier of the octree attached to 
this spatial frame. 

10. Shape contents enclosed by the outer boundary. This 
can be void when the frame defines a hole.. 

11 . Methods 

a. Put/GetNeighbor 
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b. Put/GetLogicalNeighbor 

c . Put/GetShapeContents 

d. Put/GetSurf aceContents 
e - ReseampleSurf ac.e 

5 f. Put /Ge tout erBoundary 

g . Put/GetlnnerSpatialFrameList 

h J Put /GetParentSpat ialFrame 

i . Put /GetTentporalFrame 
12 . Evolutionary law 

10 i a. Parameter list 

b. Boundary conditions 

c. Initial conditions 

d. Status 

15 Shape index manifold 

1 . Database identifier 

2 • Status 

3 . Floating point references 

20 a. XYZ 
b. ABC 

■ 4. Misfit threshold 

5-. Shape index charts 

6. Boundary representation node 



25 



Shape index chart 



1. Database identifier 
2 . Symmetry group isomorphism class 
30 {Sphere, Cylinder, Cone, Frustum, Torus, 

Ellipsoid, Hyperboloid_JType_I, 
Hyperboloid_ Type_II, Circle, Ellipse, Cube, 
- Tetrahedron, Square, Triangle, Rectangle, Cyclide, 
Trivial} 
35 3 Parameter vector ABC 

4 . Canonical position transformation [Rotation matrix 
and translation vector] 
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. 5. Shape index- manifold 

6 . Shape index interval 

7. Critical point list 

8. (Static) analytical methods 



5 


a . 


Evaluator 




b. 


N.ormal 




c . 


Tangent 




d. 


First Fundamental Form 




e. 


Chris tof f el symbols 


10 


f . 


Second Fundamental Form 




g. 


Gaussian curvature 




h. 


Mean curvature 




i . 


Principal curvature 




• 

3 . 


Shape index 


15 


k. 


Signed distance function 




1. ' 


Closest point transform 




in. 


Parameter estimation 




n. 


Extrude 




o. 


. Offset 


20 


P. 


Adap t i ve s amp 1 ing 




9 . Standard methods 




a. 


Sve 




b. 


Restore 




c. 


I/O compress 


25 


d. 


I/O decompress 




e. 


Delete 




f . 


Copy 




g. 


Newa 




10. 


Misfit reduction 


30 


a. 


Absolute amount 




b. 


Relative amount 




c . 


Base type 



{Ellipse, Circle, Square, Triangle, Rectangle} 
d. Parameter vector list 
35 {u, v, normal offset} 
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ID Boundary 

1. Conic section class 

{Ellipse, Circle, Line, Parabola, Hyperbola 
5 2 . Plus side 

a. Shape index chart 

b. UV start stop 
3 . Minus side 

a. Shape index chart 
.10 b. UV start stop 

4. Boundary representation node 

Critical point assertion 

15 (We specify two lists, separating min/max from saddle 
point . ) 

1. Shape index chart 

2 . UV location 

20 3. Chart symmetry 

We use Dustin Lister's trick of using a 2s 
complement 16-bit integer to subdivide length scale field 
data into (+/-) 32K increments. Given a reservoir model 

25 of size 9 square miles by 50000 feet, we can resolve data 
to 1 cm. Applying Dustin' s idea to Shape index charts and 
dependent lists, we obtain a surprisingly compact data 
structure. For example, suppose that a shape index 
manifold contains 48 charts such that each chart supports 

30 48 misfit reduction u boxes". Further, suppose that each 
chart contains 4 critical point orbits relative to the 
chart symmetry group. : Then the Shape Index manifold will 
consume on the order of 34000 bytes « 33 Kbytes of - 
memory. These choices are arbitrary, but far exceed the 

35 empirically estimated values for the Top of Salt surface 
in the November 2002 Implicit Surface Proof of Concept 
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5 



data set. In that demonstration, this surface consumed 
slightly more than 1 MByte of memory, which implies that 
curvature-based methods are naturally 30 times more 
conservative in memory utilization. 

Quadratic Taylor expansion disk class 



1 . Status 

{ CRITIC AL_POINT_EXPANS ION} 

10 2 . Expansion point 

3 . Principal curvatures 

4. Rotation matrix and translation vector 

5. Error tolerance 

6. Shape index interval 

15 7 . Trimming curve discrete samples list 

Characteristic length scale cone class 

1. Equivalence class 

20 2. Status 

3 . . Length scale 

4 . Top 

a . Radius 

b . Center 

25 c. Percent error 

5 . Bottom 

a . Radius 

b . Center 

c. • Percent error 
30 6. Slant height ratio 

•-. An example of applying • the techniques described 
above to a relatively complicated surface will now be 
discussed. The example used herein is a synthetic 
35 surface designed to be an w unf riendly" test case. 
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Approximately 14500 triangles were created from 7349 
vertices. The representation requires approximately 1Mb 
of memory. 

Note that discrete estimators for Gaussian and mean 
5 curvature were used to compute the shape index. See, M. 
Desbrun, MJ Meyer, P. Schroeder, A. Barr, Discrete 
Differential-Geometry Operators in nD , Technical Report 
Caltech, 2000. For example, the discrete estimator for 
mean curvature H*(P) evaluated at a point P is 

10 

H *(P) = S {cota t - cot P t )-(P-Q t ) 

4 *, area * 

Figure 42 is a schematic illustration of the 
definition of the variables in the mean curvature 
estimator. See, G. Stylianou, Compariso n of triangulated 
15 surface smoothing algorithms , 

http://3dk.asu.edu/archives/seminars/ 

presentation/Compsmthal .ppt . 

The discrete Gaussian curvature estimator K*(P) 
evaluated at a point P is given by 

K * (P) = 2tu - 2 9 t (P)> 
i 

where 0, is the i ,h triangle's angle in the star of P rooted to 1 

Figure 43 is an image of the co- triangulated 
25 irregularly spaced mesh that represents the present 
example of the top of the salt. 

Several smoothing methods for triangulated surfaces 
have been considered, but all of them assume a regular 
spaced carrier grid. See, A. Hubeli, Multireso lution 
30 Techniques for Non-Manifolds, D iss ETH No. 14625, 

Computer Graphics Laboratory, Department of Computer 
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Science, ETH Zurich, 2002. This surface is sampled 
irregularly, the' mesh surface is re-samples by 
interpolating missing sample points. Figure 44 shows the 
surface when re-sampled by interpolating missing sample 
5 points. The surface of Figure 44 has approximately 10000 . 
points versus 7900 points in the original mesh shown in 
Figure 43 . 

The surface of Figure 44 is noisy everywhere, but 
the -noise is low amplitude. The . preferred smoothing 
10' "software according to A. Hubeli, Multiresolution 
Techniques for Non-Manifolds, D iss ETH No. 14625, 
Computer Graphics Laboratory, Department of Computer 
Science, ETH Zurich, 2002, incorporated herein by 
reference. The smoothing software moves points in the 
15 horizontal plane to minimize curvature, so it is 

necessary to project one version of the surface onto 
another when different numbers of iterations are applied 

Figure 45 shows the re-sampled surface after 25 
iterations of smoothing. This surface is significantly 
20 smoother, but after comparing the two images we concluded 
that we have not sacrificed its intrinsic low frequency 
contents. The mean curvature flow preferentially removes 
high frequency noise and leaves intact low frequency 
intrinsic shape. 
25 Figure 46 shows the shape index map for the re- 

sampled triangulated surface after 25 iterations. There 
are regions of approximately the same shape index that 
contain islands that have significantly different shape 
index. "An example" of this behavior is the red body in the 
30 bottom middle that contains two blue bodies. 
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We consider the following labeling options for this 
image. We recall that "label" is a grid cell equivalence 
class identifier. 

5 1. Label each blue body and the surrounding red body as 
separate shapes . 
2 . Label the two blue bodies as part of the same shape 
and consider the blue shape to be distinct from the 
enclosing red shape. 
10 3. Label the red body as blue body misfits. 

4. Label the blue bodies as red body misfits. 

Figure 47 shows that there are 33 shapes in the 
example. According to the data, structure definition 

15 described above, about 512 bytes of storage are needed to 
express a shape with zero misfit reduction. Assuming 
adequate misfit reduction per shape consumes another 512 
bytes, we conclude that we can well approximate this 
noisy surface in 33 Kbytes of data structure. We need 

20 about 1 Mbytes of storage to define this surface, so we 
obtain a 1000 / 33 = 3 0-fold compression. 

It is important to understand that we do not equate 
a gepmetrical label to a geological theory. We prefer a 
description of the topography that uses the least number 

25 of parameters. 

According to preferred embodiments of the invention 
an important implicit surface data structure is the 
narrow band octree. The curvature analytical 
30 enhancements to implicit surface technology described" 
herein enable the construction a reservoir framework 
using significantly less memory and CPU than conventional 
triangulated surface technology requires. Additionally, 
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the curvature-based methods describe herein enable the 
construction of a new generation of structural framework. 
We enumerate a number of advantages of curvature analysis 
in support of implicit surface methods, according to 
preferred embodiments of the present invention. 

I. Computer resource economy 

a. We describe the shape of a geological surface 
in terms of explicit surface forms, e.g., 
ellipsoid and torus, rather than contextually 
in terms of a technology used to represent 
shape , e.g., triangles . 

b. We replace numerical computation of curvature, 
normals, signed distance functions, etc. by 

5 analytical formulae which enables us to drop 

the buffering needed to hold the sampled data. 

II • Surface flattening 

c. Many operations on surfaces are easier to 

0 perform when a flattened version of a surface 

is available. 

d. We refer to the flattening as the 
"parameterization" of the surface. 

e. - When we construct the parameterization of a 

15 surface, then we also generate an estimate of 

the geometric strain needed to flatten the 
surface. Strain says how comparable is the 
flattened version to the given surface.. 

f. Parameterization is much more useful if the 
0 flattened surface can be placed in the 

framework, rather than on. some planar .surface 
that is outside the ' framework . Shape analysis 
enables us to embed the parameter space surface 
in our xyz coordinate space. 

>5 
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III. 3D conformal grid generation and forensic 
reconstruction 

g. Using surface parameterisation, shape analysis 
knows how extract a 3D conformal grid from each 

5 framework sub-volume. 

h. Regular grids are convenient for simulation.. 
That does not mean that we have to permanently 
record a regular sampling in regions that show 
little or no variation. We equate this 

10 variation to an estimate of curvature. On 

demand we provide regularly spaced samples in 
low curvature regions. 

i. We call this curvature-adaptive sampling. 
Curvature adaptive sampling does not waste 

15 space on recording redundant information. 

j . Curvature-dependent adaptive grid generation is 
more economical than the non- curvature -adaptive 
form. An application can easily extract a 
perfectly uniform 3D grid on demand. . 

20 k. Mean curvature flow can be used to fill in 

holes and to construct a plausible 
representation of an eroded region of a 
sequence. 

25 IV. Framework construction, editing, and reconstruction 

1. Based on the implicit surface's signed distance 
function, we construct a connected sum of 
shapes that are sectors of a conic section 
fibre bundle. 

30 m. Boolean operations that take shape into account 

are more efficient than those that do not. 
n. Sending and receiving elementary shape 

descriptions is inherently more efficient than 
sending and receiving the corresponding 
35 triangulation. The Marching Cubes algorithm can 

be used to triangulate an implicit surface's 
. elementary shape constituents. 
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o. The. identification of elementary shapes is 
based on application of mean curvature flow 
(mcf ) . 

p. We <do not worry if a process is well 
5 approximated by mcf- Rather all that matters is 

whether the elements in a framework can be 
defined as the solution to a flow problem, 
q. We use structure synopsis diagrams to send and 
receive a high level resume of the connectivity 
10 and shape of a structural framework. This data 

structure extends the traditional Reeb graph 
summary of connectivity to a shape description 
that is intuitive - which is something a byte 
stream describing millions of triangles can 
15 never hope to achieve. 

V. Noise suppression 

r. Elementary shape identification can be applied 
to a noisy surface. ' 
20 s. Shape analysis provides a tool to distinguish 

noise suppression from shape decay. 

VI. Robust shape identification 

t. Suppose that an implicit surface is represented 
25 in a 3D grid with a tri-linear interpolation 

function. We have discovered necessary 
conditions on the interpolation function 
coefficients in order that a grid cell in the 
grid contain a critical point. It is very easy 
30 and fast to test these conditions, so we can 

quickly locate a surface's critical points in 
the grid. 

u. Group* theory -yields* a* robust classification of 
symmetry induced from the tri-linear 
,35 interpolation approximation of the implicit 

surface's signed distance function. This is 
used to identify the optimal elementary shape 
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candidate shape that fits a single grid cell of 
the 3D grid representation of the implicit 
surface. 

v. It enables local quantitative comparison of two 
5 • surface shapes. The process identifies regions 

on the measured surface where there is poor 
agreement and provides an estimate of the 
energy expended in the geometrical 
transformation needed to finish the package. 

10 . 

Figure 48 is a schematic illustration of a system 
for improved extraction of hydrocarbons from the earth, 
according a preferred embodiment of the invention. 
Seismic sources 150, 152,' an<4 154 are depicted which 

15 impart vibrations into the earth at its surface 166. The 
vibrations imparted onto the surface 166 travel through 
the earth; this is schematically depicted in Figure 7 as 
arrows 170. The vibrations reflect off of certain 
subterranean surfaces, here depicted as surface 162 and 

20 surface 164, and eventually reach and are detected by 
receivers 156, 158, and 151. 

, Each of the receivers 15 6, 158, and 151 convert the 
vibrations into electrical signals and transmit these 
signals to a central recording unit 170, usually located 

25 at the local field site. The central recording unit 170 . 
typically has data processing capability such that it can 
perform a cross- correlation with the source signal 
thereby producing a signal having the recorded vibrations 
compressed into relatively narrow wavelets or pulses . In 

30 addition, central recording "units may provide other." 
processing which may be desirable for a particular 
application. Once the central processing unit 170 
performs the correlation and other desired processing, it 
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communicates the seismic data to data processor 180 which 
is typically located at the local field site. The data 
transfer from the central recording unit 170 in Figure 7 
. is depicted as arrow 176 to a data processor 180. Data 
5 processor 180 can be used to perform processing as 

described in steps 300 to 316 in Figure 50 and 330 to 340 
in Figure 51, as is described more fully below. 

The seismic data collected from recording unit 170 
which is usually a relatively large data set. 
10 Importantly, according to the invention, the data 
processing to generate an efficient and robust 
subdivision of shapes representing the seismic data is 
performed on data processor 180. In this way a 
compressed,' stable, accurate representation of new data 
15 is transmitted to other processing centers. 

At the field site, a subset. 182 of a larger earth 
model 192 is provided to aid in the field activities. 
The after the subdivision of shapes, the Fragment earth 
model 182 can be updated with the newly acquired data for 

20 use locally. I 
Data processing center 190 is located away from the 
wellsite or field site, typically at an asset management 
location. In some cases data processing center 190 is 
physically distributed across a number of separate 
25 processing centers over a wide geographic region. Data 
processing center 190 integrates the subdivision of 
shapes representing the earth structures into existing 
earth model 192. The integration of both the geometry 
framework' and material field properties is preformed". 
30 While in some cases the integration of the new shape 

information can be done at the field site, this is not 
normally done due to a number of factors including (1) 



89 



WO 2004/057375 




PCT/GB2003/005395 



lack of full data set of earth model at the field. site, 
lack of computational facilities, and lack of sufficient 
local expertise. 

Once the earth model 192 is updated with the newly 
5 acquired information, updated earth model data from model • 
192 is preferably sent. back to the field site data 
processor 180. The information sent back to the field 
site preferably includes both (1) shape information to 
update geometry framework and material field properties 
10 of earth model fragment 182, and (2) decisions relevant 
to field site activities based in part on the updated 

earth model 192 . 

For example, according to one embodiment, based on 
new information integrated in earth model 192 at. data 
15 . processing center 190, improved predictions of fluid flow 
are made though the earth structures. Based on these 
improved predictions the rate of extraction from 
production well 114 is preferably altered using surface 
equipment 116 to optimize production rates from reservoir 
20 rock 120 while minimizing likelihood of undesirable 
events such as water breakthrough. 

According to an embodiment during the wellbore 
construction, based on new information integrated in 
earth model 192 at data processing center 190, improved 
25 predictions on the likelihood of structural failure of a 
wellbore 110 being drilled into reservoir rock 120 from 
drilling rig 112. Based on these improved predictions, 
the drilling plan used to drill the well 110 is 
preferably updated, for example to reduce the 'risk of 
30 wellbore stability problems, to steer the drilling 

activity more precisely to certain locations within the 
reservoir rock and/or avoid faults, etc, 
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According to another embodiment of the invention, 
data processing center 190 can more easily communicate 
geologic information from one earth model based on 
geometrical representation to another earth model based 

•5 on a different geometrical representation. In this way, 
the invention can be used to facilitate communication of 
earth model information between different models. 
.Similarly, according to another embodiment of the , 
invention, the techniques described her used to 

10 facilitate the aggregation of geologic . information from a 
number of geometrical representations of the earth 
structures. 

Figure 49 shows further detail of data* processor 
180, according to preferred embodiments of the invention. 
15 Data processor 180 preferably consists of one or more 
central processing units 350, main memory 352, 
communications or I/O modules 354, graphics devices 356, 
a floating point accelerator' 358, and mass storage such- 
as tapes and discs 351. 
20 Figure 50 shows steps in a method for • processing 

data used for hydrocarbon extraction from the earth, 
according to preferred embodiments of the invention. In 
step 300, sampled data representing earth structures is 
received. In step 310 one more symmetry transformation 
25 groups are identified from the sampled data. In step 
312, a set of Morse theoretical height field critical 
points are identified from the sampled data. In step 314, 
a plurality of subdivisions of shapes are generated such 
that when aggregated, the subdivisions accurately, 
30 efficiently, and irobustly represent the original earth 
structures. The generation in step 314 is based on the 
set of identified critical points and the symmetry 
transformation groups, according to the techniques 
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descried above and in further detail in Figure 51. In 
Step. 316, earth model data is processed using the 
generated subdi vision of shapes. In step 318 activity 
relating to extraction of hydrocarbons from a hydrocarbon 
5 reservoir is altered based on the processed earth model 
data. Various embodiments involving processing of earth 
model data and altering activity in steps 316 and 318 are 
described above with respect to Figure 48 . - 

Figure 51 shows further detail of steps in 
10 generating an efficient and robust subdivision of shapes, 
according to preferred embodiments of the invention. In 
step 330, each of the identified symmetry transformation 
groups is preferably associated with to a plurality of 
shape families. Each of the shape families preferably 
15 includes 'a set of predicted critical points. The shape 
families for each symmetry transformation group are 
preferably contained in a lookup table. In step 332 a 
shape family is selected from the plurality of shape 
families. The selection is preferably based on closeness 
20 of correspondence, or best fit, between the identified 
critical points from the original sampled data and the 
predicted critical points of the selected shape family. 

Note that each of the symmetry transformation groups 
is preferably a set of dif f eomorphisms that act on a 
25 topologically closed and bounded region in space-time 
such that under transformation the region occupies the 
same points in space. 

Each shape family preferably has an associated set 
of symmetry transformation group orbits, each orbit being 
3 0 associated with orbit information that specifies whether 
the orbit contains a predicted critical point and value 
of the Gaussian curvature of a point in the orbit. In 
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step 334, the orbit information from the set of symmetry 
transformation group orbits associated with the selected 
shape family is preferably applied to the original 
sampled data. In step 336 the result of step 334 yields 
5 a unique specification of a shape from the selected shape 
family- In step 338, each of the plurality of 
subdivisions of shapes is preferably generated by 
identifying a part of the uniquely specified shape that 
corresponds to the sampled data. In step 340, the 
10 identified parts are assembled, or sewn together, such 
that a representation of the earth structures is 
generated. 

The assembled subdivisions are advantageously more 
efficient and robust than conventional methods. For 

15 example the number of parameters in each subdivision 

times the number of subdivisions is substantially less 
then would be needed using a faceted representation 
method, and the plurality of subdivisions is more 
numerically stable than third order or higher 

20 representation. 

While the invention has been described in 
conjunction with the exemplary embodiments described 
- above, many equivalent modifications and variations will 
25 be apparent to those skilled in the art when given this 
disclosure. Accordingly, the exemplary embodiments of 
the invention set forth above are considered to be 
illustrative and not limiting. Various changes to the 
described embodiments may be made without departing from 
30 the spirit and scope of the invention. 
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CLAIMS 

1. A method for processing data used for 
hydrocarbon extraction from the earth comprising the 
5 steps of: 

receiving* sampled data representing earth 
structures; ✓ 

identifying one or more symmetry transformation 
groups from the sampled data; 
10 identifying a set of critical points from the 

sampled data; 

generating a plurality of subdivisions of 
shapes the subdivisions together representing the 
earth structures, the generation being, based at 
15 least in part on the set of identified critical 

points and the symmetry transformation groups; and 

processing earth model data using the generated 
subdivision of shapes. 

20 2. A method according to claim 1 wherein the 

identified symmetry transformation group is a set of 
dif f eomorphisms that act oh a topologically closed and 
* bounded region in space- time such that under 

transformation said region occupies the same points in 

25 space. 

3 . A method according to claim 1 wherein each of 
the identified symmetry transformation groups corresponds 
to a plurality of shape families. 

30 
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4. A method according to claim 3 wherein each of 
the plurality of shape families comprises a set of 
predicted critical points. 

5 . A method according to claim 4 .wherein the step 
of generating subdivisions comprises selecting a shape 
family from the plurality of shape families that 
corresponds to the identified symmetry transformation 
group, said selecting being based on closeness of 
correspondence between the identified critical points 
from the sampled data and the predicted critical points 
of the selected shape family. 

6. A method according to claim 5 wherein each 
shape family has an associated set of symmetry- 
transformation group orbits, some of the orbits being 
associated with critical points and other orbits being 
associated with distinguished' Gaussian curvature values. 

7 . A method according to claim 6 wherein each 
symmetry transformation group orbit of the selected shape 
family is associated with orbit information that 
specifies whether the orbit contains a predicted critical 
point and value of the Gaussian curvature of a point in 
the orbit. 

8. A method according to claim 7 wherein the orbit 
information from the set of symmetry transformation group 
orbits associated with the selected shape family "is 
applied to the sampled data thereby generating a unique 
specification of a shape from the selected shape family. 
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9 . A method according to claim 8 wherein each of 
the plurality of subdivisions of shapes is generated by 
identifying a part of the uniquely specified shape that 
corresponds to the sampled data. 

10. A method according to claim 9 wherein the 
identified parts of the uniquely specified shapes are 
assembled, thereby generating a representation of the 
earth structures. 

11. A method according to claim 10 wherein the 
generated representation is continuous. 

12 . A method according to claim 11 wherein the 
15 generated representation is smooth. 

13 . A method according to claim 9 wherein the 
uniquely specified shapes are specified using 
dif f erentiable functions including one or more of the 
20 following types: surfaces derived from conic sections, 

splines, general polynomials and trigonometric functions. 

14. A method according to claim 1 wherein' the 
sampled data is smoothed prior to said steps of 

25 identifying critical points and identifying one or more 
symmetry transformation groups. 

15 . A method according to claim 1 wherein the 
identified critical points are Morse theoretical Height 

30 field critical points consisting of the following three 
types: minima, maxima and saddle points. 
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16. A method according to claim 15 wherein said step 
of generating a plurality of subdivisions comprises 
applying a canonical homogeneous transform such that the 
number of parameters needed to uniquely describe a shape 

5 in the earth structure is minimized. 

17. A method according to claim 1 wherein the earth 
model data is geologic data, geophysical data, 
petrophysical data, mechanical earth model data and/or 

10 reservoir fluid flow data. 

18.. A method according to claim 1 wherein earth 
model data is processed such that earth models are 
updated, alternative versions of existing earth models 
15 are created, time-lapse earth models are generated and/or 
the earth model data is distributed to other earth models 
or other applications . 

19 . A method according to claim 1 wherein the 

20 sampled data represents sampled physical structure and 
material properties of 4 the earth structures. 

20. A method according to claim 1 wherein said step 
. of processing earth model data comprises making 

25 predictions of fluid flow though at least some of the 
earth structures and wherein the altered activity is 
altering the rate of extraction based on said 
predictions . 

30 21. A method according to claim 1 wherein said step 

of processing earth model data comprises predicting the 
likelihood of structural failure of a wellbore though at 
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least some of the earth structures and wherein the 
altered activity is altering the drilling of the wellbore 
based on the predicted likelihood of failure. . 

-22 . A method according to claim 1 wherein said step 
of processing earth model data comprises communicating 
geologic information relating to at least some of the 
earth structures between a first geometrical 
representation and a second geometrical representation of 
the earth structures . 

23. A method according to claim 1 wherein said step 
of processing earth model data comprises aggregating 
information from a plurality of geometrical 
representations of the earth structures and wherein the 
altered activity is based at least in part on the 
aggregated information. 

24. A method according to claim 1 wherein said step 
of processing earth model data comprises constructing an 
earth model to a user specified error tolerance using the 
generated subdivision of shapes. 

25. A method according to claim 1 wherein each of 
25 the plurality of subdivisions of shapes -is generated by 

identifying a part of a uniquely specified shape that 
corresponds to the sampled data. 

26. A method according to claim 1 wherein the step 
30 of generating a plurality of subdivisions comprises the 

steps of: 

analyzing curvature of the sampled data thereby 
generating a shape index field; and 
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identifying functions that fit the shape index 
field. 

27. A method according to claim 26 wherein the 
functions are dif f erentiable . 

28. A method according to claim 1 wherein the 
plurality of subdivisions are generated such that the 
number of parameters in each subdivision times the number 
of subdivisions is substantially less than would be 
needed using a faceted representation method. 

29. A method according to claim 1 wherein the 
plurality of subdivisions are generated such that they 
are more numerically stable than third order or higher 
representation . 

30. A method according to claim 1 wherein the 
sampled data is a faceted representation of the earth 

20 structures . - 

31. A method according to claim 30 wherein the 
faceted representation is a triangle mesh. 

25 32. A method according to claim 30 wherein the 

faceted representation is a grid. 

33 .* A method according to claim 1 wherein the 
sampled data is data measured with seismic acquisition 
30 equipment. 

34. A method according to claim 33 wherein said 
steps of receiving, identifying one or more symmetry 
transformation groups, identifying a set of critical 
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points and generating a plurality of subdivisions of 
shapes are preformed at or. near the location where the 
sample data is measured. 

5 35. A method according to claim 34 wherein said 

step of processing earth model data is performed in one 
or more locations remote from the location where the 
sample data is measured. 

10 36. A system for improved extraction of 

hydrocarbons from the earth comprising: 

a storage system adapted to receive and store 
sampled data representing earth structures; 

a processing system adapted to identify one or 
15 more symmetry transformation groups from the sampled 

data, identify a set of critical points from the 
sampled data, and generate a plurality of 
subdivisions of shapes the subdivisions together 
representing the earth structures, the generation 
20 being based at least in part on the set of 

identified critical points and the symmetry 
transformation groups; 

an earth model processing system adapted to 
processes earth model data using said generated 
25 subdivision of shapes; and 

an interface to output the processed earth 
model data to an operator. 

37. -A system according to claim 36 wherein the 
30 identified symmetry transformation group is a set of 

diffeomorphisms that act on a topologically closed and 
bounded region in space- time such that under 
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transformation said region occupies the same points in 
space . 

38. A system according to claim 36 wherein each of 
5 the identified symmetry transformation groups corresponds 

to a plurality of shape families, each of which comprises 
a set of predicted critical points. 

39. A system according to claim 38 wherein the 
10 subdivisions are generated such that a shape family. is 

selected from the plurality of shape families that 
corresponds to the identified symmetry transf ormatiori 
group, said selecting being based on closeness of 
correspondence between the identified critical points 
15 from the sampled data and the predicted critical points 
of the selected shape family. 

40. A system according to claim 39 wherein each 
shape family has an associated set of symmetry 

20 transformation group orbits, each orbit being associated 
with orbit information that specifies whether the orbit 
contains a predicted critical point and value of the 
Gaussian curvature of a point in the orbit, and wherein 
the orbit information from the set of symmetry 

25 transformation group orbits associated with the selected 
shape family is applied to the sampled data thereby 
generating a -unique specification of a shape from the 
selected shape family. 

30 % 41. A system according to claim 40 wherein each of 
the plurality of subdivisions of shapes is generated by 
identifying a part of the uniquely specified shape that 
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corresponds to the sampled data, and wherein the 
identified parts are assembled, thereby generating a 
representation of the earth structures. 

5 42. A system according to claim 3 6 as part of a 

system adapted to assist a decision making process 
relating to. extraction of hydrocarbons from a hydrocarbon 
reservoir modeled by the processed earth model data; 

10 43. A system according to claim 36 wherein the 

plurality of subdivisions are generated such that they 
are more numerically stable than third order or higher 
representation.. 

15 44. A system according to claim 3 6 wherein the • 

sample data are acquired from' the earth structures using 
seismic acquisition equipment, the storage system and the 
processing system are located at or near the location 
where the sample data are acquired, and the earth model 
'20 processing, system is located in one or more locations 
remote from the location where the sample data is 
acquired. 

25 
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